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ABSTRACT 

We study velocity moments of elliptical galaxies in the Coma cluster using Jeans equa- 
tions. The dark matter distribution in the cluster is modelled by a generalised formula 
based upon the results of cosmological TV-body simulations. Its inner slope (cuspy or 
flat), concentration, and mass within the virial radius are kept as free parameters, 
as well as the velocity anisotropy, assumed independent of position. We show that 
the study of line-of-sight velocity dispersion alone does not allow to constrain the 
parameters. By a joint analysis of the observed profiles of velocity dispersion and kur- 
tosis we are able to break the degeneracy between the mass distribution and velocity 
anisotropy. We determine the dark matter distribution at radial distances larger than 
3% of the virial radius and we find that the galaxy orbits are close to isotropic. Due 
to limited resolution, different inner slopes are found to be consistent with the data 
and we observe a strong degeneracy between the inner slope a and concentration c: 
the best-fitting profiles have the two parameters related with c = 19 — 9.6 a. Our 
best-fitting NFW profile has concentration c = 9, which is 50% higher than standard 
values found in cosmological simulations for objects of similar mass. The total mass 
within the virial radius of 2.9 Mpc is 1.4 x 10 15 h^ Mq (with 30% accuracy), 85% 
of which is dark. At this distance from the cluster centre, the mass-to-light ratio in 
the blue band is 351 h 70 solar units. The total mass within the virial radius leads to 
estimates of the density parameter of the Universe, assuming that clusters trace the 
mass-to-light ratio and baryonic fraction of the Universe, with Qq = 0.29 ± 0.1. 

Key words: methods: analytical - galaxies: clusters: individual: Coma - galaxies: 
kinematics and dynamics - cosmology: dark matter 



1 INTRODUCTION 

The Coma cluster of galaxies (Abell 1656) is one of the most 
extensively studied in our neighbourhood (see e.g. Biviano 
1998 and references therein). Starting with the seminal pa- 
per of Kent & Gunn (1982) significant effort went into dy- 
namical modelling of the cluster. In the early studies based 
on about 300 galaxy velocities only velocity dispersion was 
modelled and it was most often assumed that the mass fol- 
lows light and that the galaxies are on isotropic orbits. Mer- 
ritt (1987) showed that if a larger variety of models is al- 
lowed there is a strong degeneracy between the dark matter 
distribution and velocity anisotropy and many models can 
be shown to be consistent with the data. Without any prior 
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knowledge on the mass distribution even considering higher 
velocity moments would probably not be of much help. 

Recently, due to theoretical progress mainly by the 
means of A-body simulations, our knowledge on possible 
dark matter distributions within gravitationally bound ob- 
jects has improved significantly. There seems to be general 
agreement at least as to the behaviour of dark matter den- 
sity profiles at large radial distances (g oc r -3 ). Whether 
the inner dark matter density profile is g oc r -1 (as in the 
so-called universal profile advocated by Navarro, Frenk & 
White 1997, hereafter NFW) or g oc r~ 3/2 (as preferred 
by Moore et al. 1998; see also Fukushige & Makino 1997) 
or is flat (as suggested by the observed rotation curves of 
dwarf and low surface brightness galaxies, e.g. McGaugh & 
de Blok 1998) is still a matter of debate (a recent analysis 
by Jimenez, Verde & Oh (2003) of high resolution rotation 
curves of spiral galaxies shows that 2/3 of the sample can be 
accounted by NFW profiles, but 2/3 also with a flat core). 
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We therefore consider a generalised profile with different in- 
ner slopes and also allow for different dark matter concen- 
trations. We constrain this variety of dark matter density 
profiles by modelling velocity moments of galaxies. In addi- 
tion to constraints from the line-of-sight velocity dispersion 
profile, we incorporate constraints from the fourth velocity 
moment, the kurtosis. 

The NFW profile has been found consistent with the to- 
tal mass distribution inferred from the galaxy data combined 
from many clusters in the CN0C1 (van der Marel et al. 2000) 
and ENACS (Biviano et al. 2003) surveys. Although van der 
Marel et al. (2000) considered higher velocity moments they 
did not apply them rigorously to further constrain the mass 
distribution. Studies based upon X-rays, assuming that the 
hot X-ray emitting gas is in hydrostatic equilibrium in a 
spherical potential, usually lead to NFW-like cuspy centres 
(McLaughlin f999; Tamura et al. 2000; Sato et al. 2000). 
The studies based on gravitational lensing focus on the inner 
shape of the density profile. The slopes agree with the NFW 
prediction in some studies (e.g. Broadhurst et al. 2000) while 
one team finds a preference for a flat core (Tyson, Kochanski 
& dell' Antonio 1998, see also Williams, Navarro & Bartel- 
mann 1999). Note that the Coma cluster is too close for its 
mass profile to be probed through gravitational lensing. 

The amount of galaxy velocity as well as brightness and 
morphological type measurements for the members of Coma 
has increased over the last two decades making it possible to 
analyse separately the samples of elliptical and spiral galax- 
ies. While ellipticals appear to be in dynamical equilibrium 
justifying the application of Jeans formalism to study their 
velocity moments, most of the spirals are probably infalling 
onto the cluster. The modelling of the infall of spirals will 
be presented in the follow-up study. 

The paper is organised as follows. In Section 2 we de- 
scribe our data. In Section 3 we present our assumptions 
concerning the matter content of the cluster, i.e. the distri- 
butions of galaxies, gas and dark matter. Section 4 outlines 
our formalism for modelling the velocity moments of ellip- 
tical galaxies based on Jeans equations. Section 5 presents 
its application to constrain the parameters of our model, in 
particular the distribution of dark matter in the cluster. The 
discussion follows in Section 6. 



2 THE DATA 

We have searched the NED database for galaxies within 300' 
of RA=12h59m35.7s, Dec=+27°57'33" (J2000) i.e. the po- 
sition of the elliptical galaxy NGC 4874, well established as 
the centre of the Coma cluster (Kent & Gunn 1982). The 
galaxies were required to have velocities between 3000 and 
11000 km/s (given the velocity dispersions we will find be- 
low, this velocity range extends to i> 4 a). For the calcula- 
tion of the velocity moments and the subsequent study of 
kinematics, we remove from the list galaxy pairs and known 
members of pairs (as given by NED), as we wish to probe the 
global cluster potential but not its local enhancements. We 
then obtain a sample of 1068 galaxies shown in the upper 
panel of Figure 1 as points in the plane of velocity versus 
distance from the centre of the cluster. To determine the 
membership of galaxies in the cluster we proceed in a sim- 
ilar fashion as Kent & Gunn (1982). As can be seen from 
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Figure 1. Upper panel: 1068 galaxies selected from the NED 
database within 300' from NGC 4874 with heliocentric velocities 
between 3000 and 11000 km/s. Middle panel: 355 E-S0 galaxies, 
members of Coma. Lower panel: 163 spiral galaxies, members of 
Coma. The curves indicate envelopes of the cluster. 
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Figure 1, the members of the cluster are well separated from 
the foreground and background galaxies in velocity space. 
We have therefore selected the probable members of Coma 
as lying within the two curves shown in Figure 1, symmetric 
with respect to v = 7000 km/s, the value close to the mean 
velocity of the cluster. This procedure leaves us with 967 
galaxies. For the determination of the luminosity distribu- 
tion, we keep the members of galaxy pairs and remove a few 
galaxies for which no magnitude estimate is available (they 
may contribute to the calculation of velocity moments how- 
ever). We then proceed with the membership determination 
as before. 

It is generally believed that only early-type (E-SO) 
galaxies can be considered in dynamical equilibrium within 
a cluster in opposition to spirals which are believed to be in- 
falling (e.g. Tully & Shaya 1984; Huchra 1985). We therefore 
construct separate samples of E-SO and spiral galaxies. The 
morphological type of the galaxies has been determined by 
consulting the NED, SIMBAD and LEDA databases. Among 
the 967 galaxies belonging to Coma and selected for the anal- 
ysis of the velocity moments, we find 355 E-SOs, 163 spirals 
and 449 other galaxies for which the morphological classifi- 
cation is unknown or uncertain. As already noticed by Kent 
& Gunn (1982), the cluster shows a clear morphological seg- 
regation. The distributions of E-SOs and spirals in the R—v 
plane are shown in the middle and lower panel of Figure 1. 
Comparison of the two plots reveals that while E-SOs cluster 
and dominate in the central parts of Coma and are not nu- 
merous at radial distances larger than 80', spirals are more 
uniformly distributed and underrepresented in the central 
region. Similar subsamples are constructed for the analysis 
of the luminosity distribution. Note that iterating on the 
mean velocity and establishing new symmetric envelopes for 
the cluster in the observed phase space has virtually no ef- 
fect on the cluster membership. We note a possible group 
of galaxies at R < 10' and v < 4000 km/s that might have 
been missed, but contributes little to the central internal 
kinematics. 



3 THE MATTER CONTENT OF THE 
CLUSTER 

3.1 The mass in stars 

The mass contributed by the stars in galaxies is estimated as 
follows. Ideally, one would like to determine the surface lu- 
minosity distribution separately for different morphological 
types of galaxies, since transforming luminosity into mass 
requires mass-to-light ratios, which are known to vary with 
morphological type. However, since we do not know the mor- 
phological type for roughly half the galaxies in our sample, 
we would have to assume for these galaxies some mean value 
of the mass-to-light ratio. Combining then the fits for the 
three classes of galaxies: E-SOs, spirals and those of unknown 
type would produce large uncertainties, since the luminosity 
distribution for spirals turns out to be quite noisy. Therefore, 
in determining the total stellar mass distribution, we fit the 
luminosity distribution for all galaxies and then translate it 
to the mass distribution using a mean mass-to-light ratio. 

Magnitudes are transformed into luminosities assuming 
that all galaxies are at the same distance associated with the 
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Figure 2. The surface luminosity distribution data together with 
the best-fitting projected NFW profiles for all 985 galaxies (filled 
symbols and solid line) and the 366 known E-SOs among them 
(open symbols and dashed line). 

mean velocity of the Coma cluster. Adopting an heliocentric 
velocity of Coma of 6925 km/s (Struble & Rood 1999), and 
correcting for the velocity of the Sun with respect to the 
Local Group and for the Local Group infall onto Virgo, we 
obtain for Coma a Hubble flow velocity of 7093 km/s, which, 
for a Hubble constant of Ho = 70kms _1 Mpc -1 (assumed 
throughout this paper) gives a distance of 101.3 Mpc and a 
distance modulus of 35.03 (neglecting the peculiar velocity 
of the Coma cluster). 

As mentioned in the previous section, our sample of 
galaxies for the luminosity analysis is somewhat different 
from the one used for the calculation of velocity moments, 
as we no longer exclude the members of pairs. Now, we have 
985 galaxies that belong to Coma, among which 366 E-SOs 
and 167 spirals. The surface luminosity profile of all galaxies 
is then determined by placing 60 galaxies per radial bin. The 
resulting distribution is shown as filled symbols in Figure 2. 
The open symbols in the Figure show a similar result for just 
E-SO galaxies but this time with 30 galaxies per radial bin. 
This second distribution will be needed in the modelling 
of the velocity moments as the distribution of our tracer 
population. 

The data shown in Figure 2 do not indicate any presence 
of a core in the surface luminosity distributions, hence we fit 
them with projections of cuspy profiles. For both samples, 
of all galaxies as well as early-type ones, the distributions 
have a changing slope, so we fit them with a projection of the 
NFW profile (Bartelmann 1996; see also section 2.5 of Lokas 
& Mamon 2001). The fit is done by x 2 minimisation taking 
all points with equal weights (although we do not know the 
errors of the magnitude measurements we expect them to be 
similar for galaxies in each bin). The 3D luminosity density 
will then be 

* W = (r/r s )(l + r/r s ) 2 ' W 
where the normalisation constant h and the scale radius rs 
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are the two fitting parameters. For the sample of all galaxies 
we obtain h = 9.05 x 10 7 Z/Q/arcmin 3 and rs = 14'4. For 
the E-SOs we have h = 3.55 x lO 8 L /arcmin 3 and rs = 7' 05 
so this distribution turns out to be more concentrated and 
therefore steeper. It might be interesting to note that the 
number density distribution of the galaxies in our sample is 
less steep in the centre than the luminosity density distribu- 
tion and can be approximated by a profile with a core. The 
steeper distribution of luminosity is probably mainly due to 
the presence of bright ellipticals in the centre of the cluster 
(see Figure 1). 

Integrating the luminosity distribution (1) and multi- 
plying by the appropriate mass-to-light ratio in the blue 
band, T, we obtain the mass distribution associated with 
the stars in galaxies 

Mg{t) = 47rL*Tr| An r \ . ( 2 ) 

V r s r + r s J 

We adopt the mass-to-light ratio in blue band for E-SOs 
of T E = 8M /L© and for spirals T s = 3Mq/L & . For 
the mass distribution of all galaxies we take the weighted 
mean of T. Since, for galaxies with known morphological 
type, the E-SOs and spirals appear in the proportions of 
2.2 : 1, then assuming that the same proportion holds for 
the whole population, the mean mass-to-light ratio amounts 
to T G = (2.2T E + T s )/3.2 = 6.43M /L . 

3.2 The gas distribution 

To approximate the contribution of the gas to the mass in 
the cluster we make use of the X-ray surface brightness dis- 
tribution which can be well approximated by the following 
formula 

„-i -36+1/2 



S(r) = So 



1 + 



(3) 



Briel, Henry & Bohringer (1992) analysed the ROSAT ob- 
servations of Coma out to 100' in the energy band of 0.5 
to 2.4 keV and obtained the following best-fitting param- 
eters: So = 4.6 x 10~ 13 erg cm" 2 s _1 arcmin~ 2 , b = 0.75 
and the core radius r c = 10.'5. Assuming that the radiation 
is produced by bremsstrahlung this distribution is obtained 
by projecting the 3D emissivity integrated over the energy 
band and assuming that the gas is approximately isother- 
mal. Since the emissivity is proportional to the square of the 
number density of electrons in the gas n, we obtain 

2l - 36 / 2 



n(r) = n 



1 + 



(4) 



where the central electron density (with h = 0.7) is no = 
3.42 x 10" 3 cm" 3 (Briel et al. 1992; Henry & Henriksen 
1986). Integrating equation (4) and multiplying by the mass 
per electron we obtain the mass distribution associated with 
the gas 



., , , 4 , , 3 ^ / 3 36 5 

M s( r ) = -Tm (m c + jm p )r 2-F1 I -, — ; -; 



(5) 



where 2-F1 is the hypergeometric function and 7m. p with 
7 = 1.136 is the mean mass of the positively charged ion 
in the gas per unit charge (m p being the proton mass), as- 
suming the gas has the primordial composition with helium 
abundance of F He = 0.23 - 0.24. 



The formula (3) for the surface brightness is the same as 
the one appearing in the so-called /3-model (we have replaced 
(3 by b to avoid confusion with the anisotropy parameter of 
the next section) , however we do not accept all the assump- 
tions of the model here concerning e.g. the relation between 
the gas and galaxy distributions. The only assumption going 
into the derivation of the gas mass distribution in equation 
(5), besides spherical symmetry, is that the gas is isother- 
mal, which is justified by recent observations of Coma with 
XMM-Newton by Arnaud et al. (2001), who find very little 
temperature variation at least in the central region of radius 
20'. 



3.3 Dark matter 

We study different possible density distributions of dark 
matter described by the following general formula 



p(r) 



Pcha 



(6) 



(r/r s ) a (1 +r/r s ) 3 - Q ' 

where p c har is a constant characteristic density, while r s is 
the scale radius of the dark matter (in general different from 
that of the luminous distribution, rs). The profiles differ 
by the inner slope r~™ but have a common outer limit- 
ing behaviour of r~ 3 . We will consider a values limited by 
< a < 3/2, which covers a wide range of possible inner 
profiles: from very steep to core-like. The cuspy profiles of 
a > are motivated by the results of cosmological Af-body 
simulations. The profile with a = 1 corresponds to the so- 
called universal profile proposed by NFW as a fit to the 
profiles of simulated haloes, while the profile with a = 3/2 
is identical to the one following from higher resolution simu- 
lations of Moore et al. (1998). The core profile with a — is 
favoured by some observations of galaxies and clusters and 
is very similar (but not identical) to the profile proposed by 
Burkert (1995). 

The scale radius r B introduced in equation (6) marks 
the distance from the centre of the object where the profile 
has a slope equal to the average of the inner and outer slope: 
r -(3+«)/2 Qj-jjgj. parameter that controls the shape of 

the profile is the concentration 



(7) 



where r v is the virial radius, i.e. the distance from the centre 
of the halo within which the mean density is A c times the 
present critical density, p C rit,o- Although in most cosmologi- 
cal Af-body simulations A c = 200 is assumed and kept con- 
stant, the value, following from the spherical collapse model, 
depends on the underlying cosmology, e.g. A c w 178 is valid 
for the Einstein-de Sitter model, while, for the currently 
most popular ACDM model with Qm ~ 0.3 and Qa = 0.7, 
we have A c « 102 (Eke, Cole & Frenk 1996; Lokas & Hoff- 
man 2001). We will keep A c = 102 in the following. 

The concentration of simulated dark matter haloes has 
been observed to depend on the virial mass. Jing & Suto 
(2000) tested the relation c(M v ) for the masses of the order 
of normal galaxies and clusters in the case of density profiles 
with a = 1 and a = 3/2 and found concentrations slowly de- 
creasing with mass (thus confirming the original observation 
of NFW) and lower for a — 3/2 than for a — 1. The only 
study using the value of A c appropriate for a given cosmo- 
logical model is the one of Bullock et al. (2001) who found 



© 0000 RAS, MNRAS 000, 000-000 



Dark matter distribution in the Coma cluster 5 



the profiles of presently forming haloes to be well fitted by 
the NFW formula with concentrations depending on mass in 
the ACDM model with the above parameters approximately 

as 

-0.122 



c(M„) = 5.95 



M v 



W ls h- 1 M fl 



© 



(8) 



In the following, we will treat the concentration as a free 
parameter using this relation to guide us as to order of mag- 
nitude of c. 

We normalise the density profile (6) so that the mass 
within r v is equal to the so-called virial mass 

4 3 

M v = -7rr„A c p cr it,o • (9) 

The characteristic density of equation (6) then becomes 

_ (3 - a)A c p crit ,o e" 
Pchar 3F^c) ' (10) 

where F a (c) is given by the hypergeometric function 
F a {x) = 2 Fi(3-a,3-a;4-a;-x) (11) 
3 



x 



= < 
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The dark mass distribution following from (6), (9) and (10) 
is 

M DW = ^»-^, (12) 
where we introduced s = r/r v . 



4 MODELLING OF THE VELOCITY 
MOMENTS 

The purpose of this work is to constrain the distribution of 
dark matter in the Coma cluster by studying the velocity 
moments of the population of elliptical galaxies in the clus- 
ter which we believe to be in equilibrium (thus neglecting 
radial streaming motions). We defer to a later paper the 
kinematical study of spiral galaxies in the context of infall. 
We will also assume that the system is spherically symmet- 
ric and that there are no net streaming motions (e.g. no 
rotation) so that the odd velocity moments vanish. 

At second order, the two distinct moments are v 2 and 
Vg = which we will denote hereafter by a 2 and erf respec- 
tively. They can be calculated from the lowest order Jeans 
equation (e.g. Binney & Mamon 1982) 

d, 2 , 2/3 2 d$ 

_ ( Wr ) + — vaT = -v— , (13) 

where v is the 3D density distribution of the tracer pop- 
ulation and $ is the gravitational potential. We will solve 
equation (13) assuming the anisotropy parameter 



p = l- 



a 2 {r) 



(14) 



to be constant with —00 < (3 < 1. This model covers all 
interesting possibilities from radial orbits (/3 = 1) to isotropy 
(P = 0) and circular orbits (P — > —00). 

The solution of the lowest order Jeans equation with 
the boundary condition a r — > at r — > 00 for f3= const is 



va 2 (fi = const) = r 



-2/3 



f°° 213 d$ 
J " ^ 



dr 



(15) 



However, the measurable quantity is the line-of-sight veloc- 
ity dispersion obtained from the 3D velocity dispersion by 
integrating along the line of sight (Binney & Mamon 1982) 



cfos(-R) 



I(R) 



R 2 ' 



R 2 



dr , 



(16) 



where I(R) is the surface distribution of the tracer and R is 
the projected radius. Introducing equation (15) into equa- 
tion (16) and inverting the order of integration, we obtain 



Clos(-R) 



2G_ 

I(R) 



/•oc 

/ da; v 
Jr 



(x)M(x)x 



2/3-2 



(17) 



-2/3+1 



v 



R 2 



where M(x) is the mass distribution and we used new vari- 
ables x and y instead of r to avoid confusion. The calculation 
of uios can then be reduced to one-dimensional numerical 
integration of a formula involving special functions for arbi- 
trary (5 = const. 

It has been established that by studying ai oa (R) alone 
we cannot uniquely determine the properties of a stellar sys- 
tem. In fact, systems with different densities and velocity 
anisotropies can produce identical ai os (R) profiles (see e.g. 
Merrifield & Kent 1990; Merritt 1987). It is therefore inter- 
esting to consider higher-order moments of the velocity dis- 
tribution. For the fourth-order moments, the three distinct 
components v$, Vg = and v 2 Vg = v 2 v^ are related by two 
higher order Jeans equations (Merrifield & Kent 1990). 

In order to solve these equations we need additional 
information about the distribution function. We will restrict 
ourselves here to functions which can be constructed from 
the energy-dependent distribution function by multiplying it 
by a function of angular momentum f(E,L) = fo(E)L~ 2fi 
with p — const. The solution of the Jeans equation for the 
fourth-order moment (see Lokas 2002) 



d — 2p — 

— {vv$) -\ vv% 

dr r 



, q ad* 
+ 3va r — — 
dr 







(18) 



vv${P = const) = 3r~ 2/3 / r 20 vo 2 r {r)%% dr . (19) 

dr 



By projection, we obtain the line-of-sight fourth moment 
2 



vL(R) = 



where 



I(R) 



Jn Vr^ 



R 2 



g(r,R,p) dr 



(20) 



(21) 
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5 RESULTS 
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Figure 3. The line-of-sight velocity dispersion (upper panel) and 
the dimensionless line-of-sight kurtosis parameter (lower panel) 
of E-S0 galaxies. The curves represent models with stars, hot gas 
and dark matter with NFW distribution of mass M v = 10 15 Mq 
and concentration c = 6. 



Introducing equations (15) and (19) into (20) and in- 
verting the order of integration, we obtain 



< B (R) = 



6G 2 

W) JB 



-2/3+1 



R 2 



g(r,R,/3)dr 



M(p) 



(22) 



and the calculation of «j 4 os (R) can be reduced to a (rather 
long) double integral. A useful way to express the fourth 
projected moment is to scale it with af os in order to obtain 
the projected kurtosis 



Klos{R) — 



< B (R) 
< S (R) 



(23) 



We use our sample of galaxy velocities to calculate the ve- 
locity moments of the E-SO galaxies. The 355 available E-SO 
galaxies were divided into bins of 39 objects. Figure 3 shows 
the line-of-sight velocity dispersion (upper panel) and kur- 
tosis (lower panel) together with their sampling errors cal- 
culated using estimators based on Monte Carlo simulations 
(see Appendix A). We find the sampling distribution of a\ OB 
to be close to normal. In the case of kurtosis, shown in the 
lower panel of Figure 3 the values are actually calculated 
and their errors propagated from the quantity (log ftios) 1 ^ 10 
which we find to be normally distributed (see Appendix A). 
Normal sampling distributions of the estimators of both mo- 
ments and very weak correlations between them justify the 
use of standard fitting procedures of these quantities based 
on x 2 minimisation. 

Our purpose here is to reproduce the observed velocity 
moments using models described by equations (16) and (22), 
with the mass distribution given by the sum of the three 
contributions (2), (5) and (12) discussed in Section 3: 



M(r) =M G + M g +M D 



(24) 



whose value is 3 for a Gaussian distribution. 



The density profile v(r) of the tracer population of early- 
type galaxies is given by equation (1) and the surface bright- 
ness I(R) by its projection. While studying velocity disper- 
sion is useful to constrain the mass, the kurtosis is mostly 
sensitive to the velocity anisotropy. 

To give a feeling of its behaviour, we show in Figure 3 
the predictions of equation (22) (and the corresponding ones 
of equation (16)) for the three cases of radial (/3 = 1), 
isotropic (f3 = 0) and circular (/3 = — oo) orbits assuming 
dark matter distribution given by an NFW profile (equation 
(6) with a = 1) for the virial mass M v = 10 15 Mq and with 
concentration c = 6 (as suggested by formula (8) for the 
mass of this order). We see that for radial orbits the kurto- 
sis profile is convex as opposed to the concave shapes in the 
case of isotropic and circular orbits. Since our data seem to 
prefer a concave shape we can expect isotropic or tangential 
orbits to fit the data best. 

We begin by fitting the line-of-sight velocity dispersion 
profile shown in the upper panel of Figure 3. We consider dif- 
ferent values of a — 0, 1/2, 1, 3/2, and, for each of them, we 
determine the best-fitting anisotropy parameter /3 and dark 
virial mass M v as a function of concentration c. The best- 
fitting values of M v obtained are of the order of 10 Mq, 
which corresponds to dark halo virial radius r v = 88' or 2.6 
Mpc. The virialised region is supposed to lie within this ra- 
dius, so in the following analysis we discard the two outer 
radial bins. 

Figure 4 shows the results of fitting just the inner 7 data 
points in the upper panel of Figure 3 for a — 0, 1/2, 1 and 
3/2. The best-fitting dark virial mass M v (upper panel) and 
the velocity anisotropy parameter for early-type galaxies (3 
(middle panel) are shown as a function of concentration. 
The lower panel of Figure 4 shows the goodness of fit \ 2 
(we do not use X* /N, with N the number of degrees of 
freedom, because it is not obvious how many parameters 
can be estimated with this procedure). As can be seen in 
Figure 4, the biggest M v obtained is 1.5 x 10 15 Mq which 
gives r v = 100' or 2.9 Mpc. The data points are now all 
within this region and therefore there is no need for further 
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Figure 4. Results of fitting linc-of-sight velocity dispersion data. 
The best fitting parameters M v and (3 arc shown in the two upper 
panels as a function of concentration for different a. The lower 
panel gives the goodness of fit \ 2 - 





0.5 1 1.5 
dark virial mass 



10 20 30 40 
concentration 





0.5 1 1.5 

dark virial mass 



10 20 30 40 
concentration 



40 



e 30 



20 



10 





0.5 1 1.5 

dark virial mass 



0.5 1 1.5 
inner slope 



Figure 5. Cuts through the la, 2 a and 3 a probability con- 
tours in the parameter space obtained from fitting <r los and 
(log Kiog) 1 / 10 . Circles (and half-circles) indicate the best-fitting 
parameters. The mass is in units of 10 15 Mq. 



adjustments in the number of data points analysed. The 
lowest panel of Figure 4 proves that neither c nor a can be 
constrained from the analysis of velocity dispersion alone: 
X 2 flattens for large c and reaches similar values for all a for 
a wide range of c. 

As discussed in Appendix A, the sampling distributions 
of a los and (log Kios) 1 ^ 10 are independent, hence we can use 
the same \ 2 minimisation scheme to find joint constraints 
following from fitting both quantities. Using the total of 
14 data points at R < 80', we jointly fit our four param- 
eters, M v , (3, a and c. Contrary to the case when only veloc- 
ity dispersion was studied, the minimisation procedure now 
converges. The minimum is found at M v = 1.2 x 10 15 Mq 
(corresponding to dark matter virial radius r v — 92' = 2.7 
Mpc), (3 = -0.13, a = and c = 19 with X 2 /N = 6.1/10. 
For a better visualisation of the constraints obtained for our 
4 parameters, we plot in Figure 5 the cuts through the 4- 
dimensional confidence region in all six possible planes with 
probability contours corresponding to la (68%), 2a (95%) 



and 3a (99.7%) i.e. A X = X 



Xmin 



= 4.72,9.70,16.3, 



where Xmin = 6.1. 

Although the cuts do not tell everything about the con- 
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Figure 6. The best-fitting line-of-sight velocity dispersion (upper 
panel) and kurtosis (lower panel) profiles of E-S0 galaxies. The 
parameters of the models arc listed in Table 1. The data are the 
same as in Figure 3, except that only the 7 inner data points are 
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Figure 7. The best fitting dark matter profiles with different 
inner slopes (upper panel) and their effective slopes (lower panel) 
as a function of radial distance (in units of the virial radius). The 
thick solid line in each panel shows the same quantities for the 
standard NFW profile with c = 6. 



Table 1. Best-fitting parameters from joint analysis of <Ti oa and 
Ki os of E-S0 galaxies. M v is in units of 1O 15 M0. 



a 


M v 


P 


c 


x 2 





1.2 


-0.13 


19 


6.1 


1/2 


1.2 


-0.14 


14 


6.2 


1 


1.2 


-0.16 


9.4 


6.4 


3/2 


1.2 


-0.21 


4.9 


6.9 



fidence region, Figure 5 can be used to draw a number of 
interesting conclusions. The most striking is the behaviour 
of the confidence region in the c — a plane shown in the 
upper right corner. Its shape shows that there is a strong 
degeneracy between the two parameters and indeed almost 
equally good fits can be obtained for the values of the in- 



ner slope other than a = 0. The best-fitting values of the 
remaining parameters (together with the corresponding \ 2 
value) when different a are assumed are listed in Table 1. 
The results confirm that indeed a = gives the best fit, but 
other inner slopes cannot be excluded. The steeper the inner 
slope (the higher the value of a), however, the lower is the 
concentration required to provide good fits to the moments. 

The velocity moments obtained with the sets of parame- 
ters listed in Table 1 are shown in Figure 6 together with the 
data; they overlap almost exactly. The dark matter profiles 
following from equation (6) with the parameters from Ta- 
ble 1 are plotted in the upper panel of Figure 7. In the lower 
panel we also show the logarithmic slopes of the profiles. 
In both panels our best-fitting profiles are compared to the 
"standard" NFW profile with concentration c = 6 (as sug- 
gested by formula (8) for the mass of the order of 10 15 Mq). 
As can be seen in the upper panel, all our profiles overlap 
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Figure 8. Upper panel: The mass distributions for the best- 
fitting model (a = 0). Lower panel: The mass-to-light ratios for 
the best-fitting models with different a. The dark matter virial 
radii are r v = 2.7 Mpc for all models shown. 

in a wide range of radial distances 0.03 r v < r < r v , and 
are somewhat steeper beyond 0.05 r v (see the lower panel) 
than the standard NFW model deduced from cosmological 
./V-body simulations. Indeed, our best fitting NFW profile 
has concentration c = 9.4, higher by 50% than the standard 
NFW profile for which c = 6 according to equation (8). 

The remaining parameters M v and (5 are better con- 
strained. The three panels on the left in Figure 5 show the 
probability contours for M v in the planes with the other 3 
parameters. We find the best estimate for the dark mass 
to be M v = (1.2 ± 0.4) x 1O 15 M (la error-bars). The 
anisotropy parameter f3 is very close to isotropic, the best fit- 
ting value (3 — —0.13 gives oe = 1.06 oy, i.e. the best-fitting 
orbits of early-type galaxies are very weakly tangential, al- 
though fully consistent with isotropy, while radial orbits are 
excluded at the 2 a level. 

The upper panel of Figure 8 compares the mass distri- 
butions in stars, gas and dark matter with the total mass 



distribution for our best fitting model with a = 0. At the 
dark matter virial radius r v = 2.7 Mpc, the total mass is 
1.4 x 10 15 M Q . The mass in galaxies is only 2% of the to- 
tal mass, the mass in gas is 13% of the total (the galaxies 
thus represent less than one-seventh of the baryonic mass), 
and the dark matter contributes the remaining 85% of the 
total mass. Therefore, the virial radius for the total mass is 
0.85~ 1/3 = 1.06 times that of the dark matter, i.e. 97' or 2.9 
Mpc. 

The cumulative mass-to-light ratio, i.e. the ratio of the 
total mass distribution to the luminosity distribution in 
galaxies, is 

M/L B = . (25) 

where M(r) is given by equation (24) with (2), (5) and 
(12), while L(±(r) — Mg(t")/T. In the lower panel of Fig- 
ure 8 we show M/Lb for our best-fitting models for dif- 
ferent a, with parameters from Table 1. The models differ 
towards the centre and only for a = 1 does the distribution 
tend to a constant value in this limit, since we have used 
the NFW profile to fit the luminosity distribution of galax- 
ies in Section 2. At radial distances larger than 0.3 Mpc, 
the cumulative mass-to-light ratio decreases slowly to reach 
M/Lb ~ 351Mq/Lq at the total mass virial radius 2.9 
Mpc. 



6 DISCUSSION 

We studied the velocity moments of early-type galaxies in 
the Coma cluster and used them to constrain the distri- 
bution of dark matter and velocity anisotropy. Our analy- 
sis differs from previous analysis of optical data (e.g. The 
& White 1986; Merritt 1987; den Hartog & Katgert 1996; 
Carlberg et al. 1997; van der Marel et al. 2000; Biviano et 
al. 2003; Biviano & Girardi 2003), in that 

(i) we have, for a single cluster, a larger sample of galax- 
ies, which, given their early morphological type, should be 
in dynamical equilibrium in the cluster potential; 

(ii) we remove pairs from the computation of the velocity 
moments; 

(iii) we include kurtosis in the analysis; 

(iv) we model dark matter distribution using a gener- 
alised formula inspired by the results of cosmological N- 
body simulations; 

(v) we include hot gas in our analysis. 

In comparison to studies based upon stacking of many 
clusters, our analysis of the Coma cluster benefits from not 
having to introduce errors in any stacking procedure, and 
from a cleaner removal of interlopers. On the other hand, 
the analyses of stacked clusters have the advantage of aver- 
aging out particular inhomogeneities of individual clusters 
such as Coma, expected in hierarchical scenarios of struc- 
ture formation. Indeed, the Coma cluster is known to have 
irregular structure both in projected space (Fitchett & Web- 
ster 1987; Mellier et al. 1988; Briel et al. 1992) and velocity 
space (Colless & Dunn 1996; Biviano et al. 1996). In par- 
ticular, the cluster has two central cD galaxies, NGC 4874 
and NGC 4889, of which the first one is the central galaxy of 
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the main cluster and the second probably belonged to a sub- 
cluster which has recently merged with the main cluster (e.g. 
Colless & Dunn 1996). There are other subgroups, such as 
the one associated with NGC 4839 at around 40' from NGC 
4874, close enough to the cluster centre to contribute to our 
analysis. 

The question whether the E-SO sample is relaxed and 
how the existing substructure may affect our results can only 
be fully addressed by cosmological iV-body simulations in- 
cluding galaxy formation, where all 3-dimensional informa- 
tion would be available. Although such an analysis has not 
yet been performed, the effect of the incomplete virialisation 
of structures of dark matter particles seen in cosmological 
simulations on the estimates of the mass of a single cluster 
through the Jeans equation has been addressed by Tormen, 
Bouchct & White (1997). They showed (see the bottom row 
of their Fig. 17) that even for significantly perturbed haloes 
the mass M(r) at distances larger than 2% of the virial ra- 
dius inferred by the proper Jeans analysis is within 30% 
(r.m.s.) of the true mass and departs from it by less than 
20% (r.m.s.) for average or relaxed haloes. Since the dark 
matter distribution is known to possess more substructure 
than is observed in the galaxy distribution (cosmological 
simulations predict many more Milky Way satellites than 
are observed, see e.g. Moore et al. 1999) and since struc- 
tures cease to grow sooner in flat universes with a cosmo- 
logical constant, in comparison with analogous structures 
growing within the Einstein-de Sitter model (assumed in 
the simulations by Tormen et al.), they are more regular to- 
day (see Thomas et al. 1998) than shown in the first three 
columns of the bottom row of Figure 17 in Tormen et al. We 
therefore believe that the discrepancy between our derived 
mass and the true cluster mass due to substructure and de- 
partures from equilibrium is significantly smaller than the 
uncertainty due to sampling errors of the velocity moments. 

Our results for the dark and total mass of the clus- 
ter are consistent with previous estimates. Using a combi- 
nation of X-ray and optical data Hughes (1989) found for 
his preferred model a total mass within 3.6 h^ Mpc (where 
we used the notation Ho = 70 /170 km s _1 Mpc -1 ) to be 
(1.3±0.2) x 10 15 hjQ 1 Mq, but a much wider range of masses 
if more general mass distributions were allowed. Briel et al. 
(1992) using ROSAT observations of Coma and assuming 
hydrostatic equilibrium of the gas derived a mass within the 
same radial distance to be (1.3 ± 0.4) x 10 15 Mq. Also, 
by analysing the infall patterns around Coma, Geller, Diafe- 
rio & Kurtz (1999) find a mass of (1.7±0.4) x 10 15 Mq 
within the same distance to the cluster centre as above. Ex- 
trapolating our results to this radial distance, we find an 
enclosed mass of (1.6 ± 0.5) x 10 h^ Mq, in good agree- 
ment with these three earlier estimates. 

We find a strong degeneracy between the inner slope 
and the concentration of the dark matter profile, with many 
combinations of the two reproducing our velocity profiles al- 
most equally well. In the range of inner slopes < a < 3/2 
we find that the best-fitting models have the two param- 
eters related almost linearly as c = 19 — 9.6 a with very 
little variation in the remaining fitting parameters, M v and 
f5. The particular shape of the degeneracy between c and a 
is due to the specific properties of the family of dark mat- 
ter profiles used in this analysis, coupled with our lack of 
velocity data at radii smaller than the scale radius (r s ) of 



the dark matter. Using smaller radial bins would allow us to 
probe smaller radial distances, but at the expense of larger 
sampling errors in the velocity moments. As is clear from 
Figure 7, the best-fitting dark matter profiles obtained for 
different inner slopes are almost the same for a wide range 
of distances (larger than 3% of the virial radius) and this 
can only be achieved with the combination of the parame- 
ters as given in Table 1. Note that our best fit to the data 
for the flat dark matter density profile might mean that the 
data may be even better fitted with the unphysical model 
where the dark matter density profile rises with radius in 
the cluster centre, i.e. a < in equation (6). 

Our best-fitting NFW profile (a = 1) has a concentra- 
tion c = 9.4 (in agreement with the above formula), 50% 
higher than c = 6 found in iV-body simulations. It should 
be kept in mind, however, that the concentration parame- 
ters in the simulations are subject to substantial scatter and 
that the formula (8) (giving c = 6) at masses of the order of 
10 5 'Mq is actually an extrapolation of results obtained for 
lower mass haloes (see Bullock et al. 2001). However, the 
"standard" NFW model with c = 6 is still within our la 
confidence region in Figure 5. 

In comparison, fits of the NFW profile to X-ray data 
of clusters, assumed isothermal, by Ettori & Fabian (1999), 
rescaled by Wu & Xue (2000) , as well as by Sato et al. (2000) , 
both yield c ~ 4 for the mass that we find for Coma within 
the virial radius (the latter after correction from A c = 200 
to A c = 102), which is only within our 2 a confidence region. 
It may be that the assumption of isothermal hot gas causes 
a lower concentration parameter. 

In an analysis similar in spirit to ours, Biviano et al. 
(2003) report c = 4 ± 2 for their stacked ensemble of 59 
ENACS clusters, fitted to an NFW model for the total mass 
density. Similarly, Biviano & Girardi (2003) show that NFW 
models with c = 5.5 are consistent with a stacked ensemble 
of 43 2dFGRS clusters. Moreover, the ENACS and 2dFGRS 
clusters analysed in both studies are on average less massive 
(i.e. with lower velocity dispersions) than Coma, and given 
that cosmological simulations find a decreasing c for increas- 
ing mass, the discrepancy with our result is even stronger. 

The difference of these four studies with our result may 
be explained by the fact that the authors mentioned above 
fit the total mass density profile to the NFW form, while 
our fit was done for the distribution of the dark component 
only. If the gas is distributed similarly in the X-ray, ENACS 
and 2dFGRS clusters as in Coma, i.e. it has a flat inner 
core, and the gas on average contributes a substantial part 
of the total mass, fitting the total mass distribution may 
result in flatter profiles than the dark haloes really have. 
The discrepancy may be also caused by the exclusion of 
central cDs by Biviano et al. (2003), although we feel that 
removing one or two galaxies from our central bin of 39 
galaxies will not affect our results. Alternatively, their mass 
model (obtained with the assumption of isotropic orbits) 
may not be coherent with the kurtosis profile of their stacked 
cluster. 

Our best fit models have orbits that are very close to 
isotropic, as is expected, in the central regions, where the 
two-body relaxation time for the galaxy system is consider- 
ably smaller than the age of the Universe, and also because 
cosmological simulations indicate that dark matter parti- 
cles typically have isotropic orbits in the centres of clusters 
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(Thomas et al. 1998; Huss, Jain & Steinmetz 1999). In other 
words, we find no significant anisotropy bias for the galaxies 
relative to the expectations of isotropy for the dark matter. 
Note that if we force isotropic orbits for the galaxies, we 
then obtain similar constraints for the inner slope, concen- 
tration parameter and mass of the dark matter component 
within the virial radius as those shown in Figure 5. 

One might wish to better reproduce the kurtosis profile. 
It is clear from Figure 6 that, while the velocity dispersion 
profile is well reproduced, there is room for improvement for 
the kurtosis profile. However, as discussed in Section 4, the 
kurtosis is mainly sensitive to the velocity anisotropy, which 
was modelled here by a single constant parameter. Besides, 
the curves shown in Figure 6 are the results of joint fitting of 
both moments and do not aim at reproducing the kurtosis 
alone. In spite of this, the inclusion of the kurtosis in our 
analysis allowed us to constrain the velocity anisotropy and 
other parameters of the model. 

From our best estimate of the mass distribution, the 
baryons (galaxies and gas) contribute 15% of the total mass 
at the dark matter virial radius r v = 2.7 Mpc (see Figure 8). 
If we assume that the cluster content is representative of the 
Universe as a whole, we can use this baryon fraction to esti- 
mate the density parameter fio (see e.g. White et al. 1993). 
Taking the baryonic density parameter at its currently best 
value from nucleosynthesis fi(, = 0.04 (with h — 0.7) we ob- 
tain fio = 0.26 ± 0.09 where the error comes from our 30% 
uncertainty in the dark virial mass value, 20% uncertainty 
in the gas mass (as estimated by White et al. 1993) and 10% 
error in the fi 6 value (Buries, Nollett & Turner 2001). 

Similarly, we may assume that clusters are good tracers, 
within their virial radius, of the ratio of mass to blue lumi- 
nosity. Given that the closure mass-to-light ratio (critical 
density over luminosity density) in the blue band is roughly 
1100^70 to 10% accuracy, consistent with the recent esti- 
mates of luminosity density from the 2dFGRS (Norberg et 
al. 2002) and SDSS (Blanton et al. 2001) galaxy surveys 
(after correction to the blue band), our mass-to-light ratio 
within the virial radius of 351 /i7o with 30% accuracy yields 
fio = 0.32 ±0.1. 

Combining these two estimates of the density parame- 
ter, we arrive at fio = 0.29 ±0.1 in excellent agreement with 
other determinations, for example the recent value obtained 
from the WMAP CMB experiment (Spergel et al. 2003). 
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APPENDIX A: THE SAMPLING 
DISTRIBUTIONS 

The validity of modelling presented in this paper rests on 
proper estimation of the velocity moments and their errors 
from observations. If their sampling distributions tend to 
normality, the statistics obtained from a sample (e.g. mo- 
ments) can be characterised by an expectation value and a 
'standard' error (Stuart & Ord 1994). For simplest statistics, 
like the second central moment, those can be calculated ex- 
actly from the population moments assuming the population 
properties. Since we are interested here also in quantities 
whose standard errors are not easily calculated analytically, 
we resort to Monte Carlo methods. 

The most natural estimators of the variance and kur- 
tosis from a sample of n line-of-sight velocity measurements 



n 



and 



K 



(S 2 ) 2 



where 



1 n 

v = - Vi 



(Al) 



(A2) 



(A3) 



is the mean of galaxy velocities in the sample. 

To investigate the distribution of these estimators for 
our binning of galaxies, i.e. when n = 39, we ran Monte 
Carlo simulations by selecting Af = 10 4 times n — 39 num- 
bers from a Gaussian distribution with zero mean and dis- 
persion of unity. Since velocity distributions of gravitation- 
ally bound objects in general do not dramatically depart 
from a Gaussian, this is a sufficient approximation for con- 
structing unbiased estimators of moments. 

For each of the N samples we compute the statistic 6* , 
namely S or K according to prescriptions given by equa- 
tions (Al) and (A2). The Monte Carlo estimate of our statis- 
tic is then the mean of all values obtained 
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Figure Al. Example of sampling distribotions of s and k ob- 
tained from a Monte Carlo simolation of sampling from a Gaos- 
sian parent distribution with n = 39 and N = 10 4 . 



1 M 

2 = 1 

and its variance is 



var(6>*) 



N 



(A4) 



(A5) 



3 = 1 



We find that the best estimates obtained in this way are 
biased, especially for kurtosis (it is interesting to note, that 
the value of kurtosis is underestimated by a few percent even 
for much larger samples with n of the order of few hundred). 
In addition, while the sampling distribution of velocity dis- 
persion is Gaussian to a very good approximation, the one 
for kurtosis is strongly skewed. Using this knowledge we con- 
struct unbiased and Gaussian-distributed estimators of line- 
of-sight velocity dispersion s and kurtosis-like variable k 

. I/ 2 



k = 



log 



2.75 



-K 



1/10 



(A6) 



(A7) 



The factor n — 1 in equation (A6) is the well known cor- 
rection for bias when estimating the sample variance, valid 
independently of the underlying distribution. In (A7) the 
factor 3/2.75 corrects for the bias in the kurtosis estimate, 
i.e. unbiased estimate of kurtosis is K' — 3K/2.75, while 
the rather complicated function of K' assures that the sam- 
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k 



0.95 - 



1.00 - 




distributed to a very good approximation and very weakly 
correlated with \g\ < 0.07. The sampling distributions of 
s and k as well as the joint distribution look very similar 
to the purely Gaussian case. The correlation coefficient can 
increase significantly only in the presence of additional out- 
liers (more numerous than predicted by the Gaussian dis- 
tribution) in the range (2a, 3a) and (— 3a, — 2a). We have 
checked that the number of galaxies with velocities in this 
range is 1 or 2 in each bin, in excellent agreement with the 
Gaussian prediction. 



0.90 - 



We can therefore assume that, to a good approxima- 
tion, all our data points measuring velocity dispersion and 
kurtosis are independent, which justifies the use of standard 
X minimisation to fit the models to the data. 
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Figure A2. The joint distribution of s and k obtained from a 
Monte Carlo simulation of sampling from a Gaussian parent dis- 
tribution with n = 39 and Af = 10 4 . Only 2000 points arc shown 
to reduce the size of the preprint. 



pling distribution of k is approximately Gaussian. Examples 
of sampling distributions of s and k from our Monte Carlo 
simulation are shown in Figure Al. We find that the stan- 
dard errors in the case of s are of the order of 11% (in agree- 
ment with an analytic result derivable with the assumption 
of Gaussian velocity distribution) while in the case of k are 
approximately 2%. 

The measured values of <7i OB and ki os calculated from 
our velocity data using equations (A6)-(A7) and (A1)-(A2) 
are shown in Figure 3. The la error bars for the velocity dis- 
persion are 0.11s. The values of kurtosis are K' = 3A"/2.75 
with approximate la error bars propagated from the 2% 
error in k. 

It is also important to check whether the sampling dis- 
tributions of the two statistics are independent. In general, 
the covariance between, e.g. even moments derived from the 
same sample does not vanish, even for Gaussian distribu- 
tions (Stuart &: Ord 1994). However, the lowest order term 
is expected to decrease with the size of the sample n. To 
check whether n — 39 in our case is indeed large enough to 
assure independence of s and k, we construct the joint sam- 
pling distribution of the two statistics. The joint distribution 
in the form of N points with coordinates given by (s, k) pairs 
calculated from each sample is presented in Figure A2. As 
indicated by the Figure the variables are very weakly corre- 
lated, we find the correlation coefficient of \g\ < 0.02. 

To check the behaviour of s and k for velocity distri- 
butions departing from gaussianity, we repeated the Monte 
Carlo simulation, again sampling n — 39 numbers from a 
Gaussian distribution, but modifying each sample by re- 
moving the 6 most inner points and adding 3 uniformly dis- 
tributed in the range (la, 2a) and 3 in the range (—2a, —la). 
Such distributions have unbiased kurtosis estimates of the 
order of K ' — 2.2, close to the lowest value obtained from our 
data (and most strongly departing from the Gaussian value 
of 3). We find that estimators s and k are again Gaussian- 
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